library(sf)
library(tidyverse)
library(vtable)
library(texreg)
library(estimatr)
library(ggmap)
library(scales)
library(extrafont)

extrafont::loadfonts()

theme_set(theme_minimal())
theme_update(text = element_text(size = 12, family = "LM Roman 10"), 
             axis.text = element_text(size = 12, family = "LM Roman 10"))
`%!in%` = Negate(`%in%`)

setwd("/Users/hankinson/Dropbox/jpipe/replication/")

# load data
df_shp <- st_read("data_sf.shp") # add names due to column name character limits for sf files
names(df_shp)[1:33] <- c("cbsa_name", "dist_band", "GEOID", "GISJOIN", "STATEFP", "PLACEFP", "NAME", 
                         "place_name", "ideo_city", "mrp_ideology", "demshare_pres", "state_name",
                         "mentions_council", "mentions_planning", "count_council", "count_planning", 
                         "binary_council", "binary_planning", "binary_all", "count_all", "mentions_all", 
                         "prop_all", "prop_council", "prop_planning", "ztax_sf", "ztax_qa", "ztax_qa_0015",
                         "ztax_qa_1530", "ztax_qa_30xx", "latitude", "longitude", "dist_cent", "dist_band_ztax")

#### summary statistics
# generate map
usa <- st_as_sf(maps::map("state", fill=TRUE, plot =FALSE))
usa <- st_transform(usa, crs = st_crs(4326))

ggplot(usa) +
  geom_sf(color = "#2b2b2b", fill = "white", size=0.125) +
  geom_sf(data=df_shp, inherit.aes=F, shape = 1, colour = "red", fill = NA) +
  coord_sf(crs = st_crs("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"), datum = NA) +
  ggthemes::theme_map()
ggsave("figures/dist_band_ztax_map.pdf", device = cairo_pdf, width = 7, height = 6, units = "in")
dev.off()

### list of covered cities grouped by cbsa
city_list <- df_shp %>% 
  filter(count_all >= 10 & !is.na(dist_band_ztax)) %>%
  group_by(cbsa_name) %>% 
  summarise(n = n()) %>%
  arrange(desc(n)) %>%
  st_drop_geometry()

knitr::kable(city_list, "latex", col.names = c("CBSA", "Number of Places"), 
             caption = "Number of cities per CBSA in LocalView data \\label{tab:cbsa_place}") %>%
  save_kable("tables/places.tex")

# generate histograms of outcomes
# count of mentions among those with mentions
ggplot(df_shp, aes(x = mentions_all, y = ..density..)) + 
  geom_histogram(aes(y=..count..)) + xlab("Number of Mentions") + ylab ("Cities") +
  theme(legend.position="right", text = element_text(size = 12, family = "LM Roman 10"),
        strip.text = element_text(size = 12, family = "LM Roman 10"), axis.text = element_text(size = 12, family = "LM Roman 10")) 
ggsave("figures/mentions_all.pdf", device = cairo_pdf, width = 7, height = 3, units = "in")

# count of minutes
ggplot(df_shp, aes(x = count_all, y = ..density..)) +
  geom_histogram(aes(y=..count..)) + xlab("Number of Minutes") + ylab ("Cities") +
  theme(legend.position="right", text = element_text(size = 12, family = "LM Roman 10"),
        strip.text = element_text(size = 12, family = "LM Roman 10"), axis.text = element_text(size = 12, family = "LM Roman 10")) 
ggsave("figures/count_all.pdf", device = cairo_pdf, width = 7, height = 3, units = "in")

# count of mentions among those with mentions
ggplot(df_shp, aes(x = dist_band_ztax, y = ..density..)) + 
  geom_histogram(aes(y=..count..)) + xlab("Zoning Tax ($/quarter-acre)") + ylab ("Density") +
  theme(legend.position="right", text = element_text(size = 12, family = "LM Roman 10"),
        strip.text = element_text(size = 12, family = "LM Roman 10"), axis.text = element_text(size = 12, family = "LM Roman 10")) +
   scale_x_continuous(labels = label_comma())
ggsave("figures/dist_band_ztax.pdf", device = cairo_pdf, width = 7, height = 3, units = "in")

# summary statistics
labs <- c("Minutes (All)", "Binary (All)",
          "Minutes (City Council)", "Binary (City Council)",
          "Minutes (Planning Commission)", "Binary (Planning Commission)",
          "Zoning Tax (quarter area)", "Median Home Value (2020)", "Wharton Index (2018)",
          "Democratic Ideology (2016)", "Dem. Pres. Vote Share (2016)")
st(df_shp, 
   vars = c("count_all", "binary_all",
            "count_council", "binary_council",
            "count_planning", "binary_planning",
            "dist_band_ztax", "med_home_val", "WRLURI18",
            "mrp_ideology", "demshare_pres"),
   labels = labs,
   out = "latex", file = "tables/summary_stats.tex")

### bivariate relationship
tags <- c("Q1", "Q2", "Q3", "Q4")
vgroup <- as_tibble(df_shp) %>% 
  dplyr::mutate(tag = dplyr::case_when(
    dist_band_ztax < 5868 ~ tags[1],
    dist_band_ztax >= 5868 & dist_band_ztax < 26851 ~ tags[2],
    dist_band_ztax >= 26851 & dist_band_ztax < 134437 ~ tags[3],
    dist_band_ztax >= 134437 & dist_band_ztax <= 533703 ~ tags[4]
  ))
vgroup$tag <- factor(vgroup$tag,
                     levels = tags,
                     ordered = FALSE)

vgroup_count <- vgroup %>%
  group_by(tag) %>%
  dplyr::summarize(n = n(), # number in group
                   mention_mean = mean(binary_all)) # group mean

vgroup <- merge(vgroup, vgroup_count, by = "tag", all.x = T)

vgroup$binary_min <-  vgroup$mention_mean + (-qnorm(1-.05/2))*sqrt((1/vgroup$n)*vgroup$mention_mean*(1-vgroup$mention_mean))
vgroup$binary_max <-  vgroup$mention_mean + (qnorm(1-.05/2))*sqrt((1/vgroup$n)*vgroup$mention_mean*(1-vgroup$mention_mean))

ggplot(data = vgroup, mapping = aes(x = tag, y = binary_all)) + 
  geom_jitter(aes(color='blue'), alpha=0.2, height = .05) +
  stat_summary(fun.y=mean, geom="point", shape=16, size=5, color="black", fill="black") +
  geom_errorbar(data = subset(vgroup, !is.na(dist_band_ztax)), aes(  ymin = binary_min, ymax = binary_max), width=.3, position = position_dodge(.5)) +
  labs(x='Zoning Tax (Quarter-Acre) by Quartile',
       y = "Union Mention (Binary)") +
  guides(color=FALSE) +
  theme_minimal() + theme(text = element_text(size = 12, family = "LM Roman 10"),
                          strip.text = element_text(size = 12, family = "LM Roman 10"), 
                          axis.text = element_text(size = 12, family = "LM Roman 10"))
ggsave("figures/mention_on_ztax.pdf", device = cairo_pdf, width = 5, height = 4, units = "in")


### results ####
# bivariate models
lm_bi_all_dist_ztax <- lm_robust(binary_all ~ scale(dist_band_ztax), data = df_shp)
lm_bi_coun_dist_ztax <- lm_robust(binary_council ~ scale(dist_band_ztax) , data = df_shp)
lm_bi_plan_dist_ztax <- lm_robust(binary_planning ~ scale(dist_band_ztax) , data = df_shp)

# with ideology
lm_bi_all_dist_ztax_ideo <- lm_robust(binary_all ~ scale(dist_band_ztax) + mrp_ideology, data = df_shp)
lm_bi_coun_dist_ztax_ideo <- lm_robust(binary_council ~ scale(dist_band_ztax) + mrp_ideology, data = df_shp)
lm_bi_plan_dist_ztax_ideo <- lm_robust(binary_planning ~ scale(dist_band_ztax) + mrp_ideology, data = df_shp)

# with pres vote share
lm_bi_all_dist_ztax_pres <- lm_robust(binary_all ~ scale(dist_band_ztax) + demshare_pres, data = df_shp)
lm_bi_coun_dist_ztax_pres <- lm_robust(binary_council ~ scale(dist_band_ztax) + demshare_pres, data = df_shp)
lm_bi_plan_dist_ztax_pres <- lm_robust(binary_planning ~ scale(dist_band_ztax) + demshare_pres, data = df_shp)

screenreg(list(lm_bi_all_dist_ztax, lm_bi_all_dist_ztax_ideo, lm_bi_all_dist_ztax_pres), 
          digits = 3, include.ci=F)

screenreg(list(lm_bi_coun_dist_ztax, lm_bi_coun_dist_ztax_ideo, lm_bi_coun_dist_ztax_pres), 
          digits = 3, include.ci=F)

screenreg(list(lm_bi_plan_dist_ztax, lm_bi_plan_dist_ztax_ideo, lm_bi_plan_dist_ztax_pres), 
          digits = 3, include.ci=F)

# all meetings
texreg(list(lm_bi_all_dist_ztax, lm_bi_all_dist_ztax_ideo, lm_bi_all_dist_ztax_pres), 
       file = "tables/lm_bi_all_dist_ztax.tex", label = "tab:lm_bi_all_dist_ztax",
       digits = 3, include.ci=F, float.pos = "!htbp", dcolumn = T, use.packages = F,
       caption = "Effect of zoning tax on probabilty of having a union mention, all meetings.", 
       custom.coef.map = list("scale(dist_band_ztax)" = "Zoning Tax (std.)",
                              "count_all" = "Total Meetings",
                              "mrp_ideology" = "Dem. Ideology",
                              "demshare_pres" = "Dem. Pres. Vote (2016)",
                              "(Intercept)" = "Intercept"))

# council only
texreg(list(lm_bi_coun_dist_ztax, lm_bi_coun_dist_ztax_ideo, lm_bi_coun_dist_ztax_pres), 
       file = "tables/lm_bi_coun_dist_ztax.tex", label = "tab:lm_bi_coun_dist_ztax",
       digits = 3, include.ci=F, float.pos = "!htbp", dcolumn = T, use.packages = F,
       caption = "Effect of zoning tax on probabilty of having a union mention, city council meetings only.", 
       custom.coef.map = list("scale(dist_band_ztax)" = "Zoning Tax (std.)",
                              "mrp_ideology" = "Dem. Ideology",
                              "demshare_pres" = "Dem. Pres. Vote (2016)",
                              "(Intercept)" = "Intercept"))

# planning boards only
texreg(list(lm_bi_plan_dist_ztax, lm_bi_plan_dist_ztax_ideo, lm_bi_plan_dist_ztax_pres), 
       file = "tables/lm_bi_plan_dist_ztax.tex", label = "tab:lm_bi_plan_dist_ztax",
       digits = 3, include.ci=F, float.pos = "!htbp", dcolumn = T, use.packages = F,
       caption = "Effect of zoning tax on probabilty of having a union mention, city planning meetings only.", 
       custom.coef.map = list("scale(dist_band_ztax)" = "Zoning Tax (std.)",
                              "count_planning" = "Total Meetings",
                              "mrp_ideology" = "Dem. Ideology",
                              "demshare_pres" = "Dem. Pres. Vote (2016)",
                              "(Intercept)" = "Intercept"))


